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ABSTRACT 

In  regions  that  do  not  contain  good  station  coverage,  a  reasonably  well  known 
velocity  structure  and  moderate  or  large  magnitude  events  (the  latter  needed  for  calibration 
using  teleseismic  constraints  on  the  location  of  a  master  event),  it  is  very  difficult  to  obtain 
accurate  regional  event  locations.  In  this  situation  standard  error  ellipses  generated  during 
event  locations  may  give  the  impression  that  an  event  location  is  better  determined  than  it 
really  is.  False  interpretations  of  the  data  may  be  due  to  a  three  dimensional  velocity 
structure  bias  unaccounted  for  by  the  method,  or  simply  to  phase  misidentification. 

Forward  waveform  modeling  does  help  constrain  the  structure,  however  this  approach 
contains  a  number  of  pitfalls.  Modeling  works  best  if  ground  truth  is  available  at  several 
epicentral  distances  and  azimuths.  Otherwise  an  in^propriate  model  might  be  constructed 
by  matching  synthetic  waveforms  and  arrival  times  to  observed  data  assumed  to  be  at  the 
wrong  location.  That  model  might  thereafter  be  used  to  locate  additional  events.  Regional 
arrays,  or  arrays  originally  constructed  for  teleseismic  monitoring,  can  sometimes  be  used 
to  identify  regional  phases  using  frequency-wavenumber  (F-K)  measurements  of  phase 
velocity  in  short  (3-5  second)  windows  on  the  seismograms,  however  this  only  gives  a 
solution  if  the  signal  to  noise  ratio  is  high.  The  F-K  method  is  sensitive  to  the  chosen 
center  frequency  and  bandwidth,  and  to  three  dimensional  heterogeneities  surrounding  the 
array.  Some  of  these  practical  issues  are  discussed  in  scenarios  involving  events  located 
within  500  km  of  a  recording  station  or  array.  We  model  seismograms  using  reflectivity 
and  linear  finite  difference  techniques,  assuming  flat  isotropic  layers  for  simplicity. 

Phases  that  can  be  identified  and  modeled  at  these  distances  include  P  and  S  wave 
reflections  and  refractions.  The  Lg  phase  cannot  be  used  as  a  phase  with  fixed  group 
velocity  throughout  this  distance  range.  This  study  highlights  the  importance  of  developing 
regional  velocity  structures  for  the  crust  and  upper  mantle  that  are  weU  constrained  by 
controlled  source  experiments. 

Key  words:  phase  identification,  event  location,  crustal  velocity  structure,  waveform 
modeling 
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OBJECTIVES: 


Our  objectives  are  to  (1)  provide  practical  solutions  to  the  problem  of  seismic  phase 
identification  in  a  region  sparsely  populated  by  seismic  stations,  and  (2)  use  these  phases  to 
locate  seismic  events  at  regional  distances.  The  underlying  scientific  problem  is  how  to 
determine  an  appropriate  velocity  structure  for  a  region  in  the  absence  of  ground  truth  or 
large  well  constrained  sources  that  have  been  locat^  teleseismically. 

RESEARCH  ACCOMPLISHED: 

We  consider  a  generic  region  in  which  small  infrequent  events  occur.  These  might 
typically  be  earthquakes  or  explosions  related  to  mining  activity.  We  assume  that  only  one 
station  or  array  has  recorded  these  events  due  to  their  size.  That  station  or  array  is  situated 
within  500  km  of  the  sources  under  study,  therefore  the  crust  and  upper  mantle  structure 
have  a  strong  effect  on  wave  propagation.  The  question  in  a  CTBT  might  then  be  to 
determine  where  a  specific  event  X  in  that  region  is  really  located. 

The  problem  is  that  we  do  not  know  tire  velocity  structure  in  this  region.  Because 
of  the  poor  station  coverage  we  cannot  perform  a  joint  inversion  to  reliably  obtain  velocity 
structure  and  event  locations,  as  has  been  done  by  Crosson  (1976)  or  by  Goins  et  al. 
(1981). 

If  the  locations  and  origin  times  of  the  events  that  have  been  recorded  at  the  station 
or  array  are  only  approximate,  we  carmot  determine  a  velocity  structure  by  simply  fitting 
observed  phase  arrival  times  with  theoretical  travel  time  curves  as  one  would  do  in  a 
controlled  source  refraction  experiment  Fitting  S-P  times  is  not  practical  either,  because 
seismograms  plotted  on  a  travel  time  (t)  versus  distance  (x)  curve  are  free  to  move  in  both 
the  t  and  x  directions  due  to  the  unconstrained  hypocenter  and  origin  time.  If  the 
hypocenter  of  one  of  the  events  is  known  independently,  but  the  origin  time  is  not  precisely 
known,  that  event  can  be  used  to  help  construct  a  travel  time  curve.  In  that  case  the  S-P 
time  provides  a  useful  modeling  constraint  because  the  position  of  the  seismogram  along 
the  distance  axis  is  fixed. 

It  is  possible  to  assume  a  starting  model  for  the  crust  and  upper  mantle  in  the  region 
by  borrowing  a  model  from  another  better  studied  region  with  similar  tectonic 
characteristics.  Continental  crust  and  upper  mantle  structures  have  been  obtained  in  a 
variety  of  regions  around  the  world  (Christensen  and  Mooney,  1995),  and  it  is  often 
possible  to  find  a  reasonable  starting  structure  by  analogy  witii  another  region. 
Alternatively,  one  can  determine  a  starting  model  using  surface  waves  that  have  crossed  the 
region  under  study.  Once  a  starting  model  has  been  determined,  the  whole  suite  of  regional 
events  can  be  relocated  based  on  that  model.  Geologically  reasonable  end-member  models 
can  be  tested  to  search  for  some  measure  of  confidence  that  includes  the  effect  of  the 
velocity  structure.  If  we  perform  a  joint  (relative)  relocation  of  all  events  using  multiple  P 
and  S  arrivals  ,  the  smallest  RMS  error  should  correspond  to  the  best  of  the  velocity 
models  considered. 

Figure  1  is  a  set  of  regional  events  that  have  occurred  in  some  arbitrary  area.  These 
events  can  be  plotted  in  record  section  format  only  because  some  crustal  structure  has  been 
selected  and  the  events  have  been  located  accordingly.  The  lineup  of  the  P  and  S  arrivals 
gives  the  impression  that  the  structure  can  be  extracted  from  this  plot  by  superimposing 
travel  time  curves  on  the  plot  However  the  lineup  of  the  phases  has  been  artificially 
imposed  by  the  assumption  of  a  specific  velocity  structure.  Other  structures  and  resulting 
locations  wiU  give  new  positions  for  the  seismograms  on  the  plot  and  result  in  phase 
lineups  as  well.  Ideally  some  independent  information  (“ground  truth”)  will  help  in  the 
selection  of  the  best  velocity  structure. 

In  Figure  2  we  have  plotted  synthetic  seismograms  and  travel  time  curves  for  two 
models.  The  travel  time  curves  constructed  here  by  ray  tracing  assume  a  spherically 
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symmetric  homogeneotis  earth.  The  one  dimensional  ray  tracing  code  handles  gradients, 
and  is  based  on  Bullen’s  (1976)  treatment  of  seismic  rays  in  a  spherically  stratified  earth 
model.  We  have  chosen  a  surface  source,  and  include  direct  rays,  refracted  rays  and 
simple  reflections  off  the  Conrad  discontinuity  and  the  Moho.  We  have  not  plotted  phase 
conversions  such  as  PS  or  PmS  in  the  travel  time  curves;  they  are  included  in  the  synthetic 
seismograms.  Model  A  is  composed  of  the  upper  few  hundred  kilometers  of  the  IASP91 
earth  model  (Kennett,  1991;  Kennett  and  Engdahl,  1991).  Kennett  discussed  the 
limitations  of  this  model  recently  (Kennett,  1995),  pointing  out  that  IASP91  is  a 
convenient  hybrid  model,  but  that  it  was  constructed  in  part  by  incorporating  some  artificial 
constraints.  Furthermore  it  is  based  on  data  from  a  number  of  different  regions,  and  may 
not  be  appropriate  for  the  region  under  study.  The  IASP91  model  is  characterized  by  a 
simple  two  layer  crust  over  an  upper  mande  with  a  very  weak  positive  gradient.  Model  B 
is  generated  by  replacing  the  2  layer  crust  in  the  IASP91  model  by  a  4  layer  crust  that  is 
more  like  a  crust  in  a  continental  shield  area.  The  additional  layers  add  to  the  complexity  of 
the  travel  time  curve.  Another  model  to  test  might  be  a  continental  shield-type  crust 
overlying  a  corresponding  shield-type  mande.  Tlie  travel  time  curves  only  provide  expected 
arrivd  times  for  comparison  with  the  data. 

To  incorporate  amplitude  information  we  turn  to  waveform  modeling.  The  velocity 
structure  can  have  a  strong  effect  on  amplitudes.  For  example  a  strong  velocity  gradient  at 
the  Moho  can  boost  up  the  amplitudes  of  the  Pn  and  Sn  phases.  An  appropriate  value  for  Q 
must  be  determined  as  well.  Qp  and  Qs  can  be  determined  by  trial  and  error,  or  by 
independent  spectral  or  coda  decay  studies.  For  each  model  we  construct  synthetic 
seismograms  using  the  linear  finite  difference  and  reflectivity  methods.  The  reflectivity 
code  was  provided  by  Harley  Benz  of  the  U.S.  Geological  Survey  (Fuchs  and  Muller, 
1971;  Benz  et  al.,  1990).  The  linear  finite  difference  code  we  use  was  originally  developed 
at  Teledyne  Geotech’s  Alexandria  Laboratory  (cf.  Jih,  1993),  and  has  been  modified  to 
include  a  prototype  attenuation  operator  currentiy  being  test^  at  Phillips  Laboratory.  In 
both  cases  we  use  an  explosion  source  at  the  surface.  After  performing  the  computation  of 
displacement  seismograms  we  differentiate  to  get  velocity  records  and  apply  a  high  pass 
filter  to  simulate  short  period  seismograms.  All  calculations  are  done  on  a  Sun 
Microsystems  SPARC20  workstation.  Finite  difference  results  are  shown  in  Figure  2.  The 
finite  difference  calculation  is  computationally  intensive  and  currently  only  practical  out  to 
250  km  distance  for  the  frequency  band  of  interest.  In  Figure  3  we  have  carried  the 
reflectivity  calculation  out  to  5(X)  km  for  the  second  model.  To  avoid  instabilities  and  speed 
up  the  calculation  we  have  provided  an  upper  frequency  limit  of  2  Hz  in  Figure  5.  This 
simplifies  the  appearance  of  the  seismograms  to  some  extent  Even  with  this  frequency 
limitation  it  is  useful  to  compute  seismograms  out  to  500  km  because  the  350-500  km 
distance  range  helps  us  constrain  the  velocity  structure  of  the  upper  mantle. 

If  array  data  are  available  we  can  determine  apparent  velocities  of  regional  phases 
using  frequency-wavenumber  (F-K)  spectral  analysis.  Unfortunately  it  appears  that  there 
are  some  problems  inherent  in  this  approach,  most  likely  related  to  heterogeneity  of  the 
crust  and  resulting  frequency  dependence  of  scattering.  An  example  is  shown  in  Figures 
4a-4b.  This  event  is  a  quarry  blast  from  Blasj0  Quarry,  recorded  at  the  NORESS  array 
(see  Dainty  and  Toksoz,  1990  for  more  information),  fa  this  case  the  distance  and  azimuth 
of  the  quarry  are  known  (324  km,  243°;  location  from  Dysart  and  Pulli,  1987).  We  have 
selected  two  subsets  of  array  stations.  Rather  than  making  use  of  the  full  complement  of 
NORESS  stations  we  are  d^imating  the  array,  to  see  if  an  array  with  spacing  more 
appropriate  for  teleseismic  moititoring  can  effectively  see  the  seismic  wavefield 
decomposed  into  plane  waves  at  regional  distances.  We  seem  to  get  strong  apparently 
stable  spectral  peaks,  but  the  phase  velocities  measured  at  two  different  frequencies  are 
quite  different  from  each  other.  Neither  set  approximates  values  expected  for  a  reasonable 
flat  layered  structure.  We  start  by  selecting  6  array  stations,  with  approximately  2  km 
spacing  between  stations,  and  a  center  frequency  of  1  Hz  and  bandwith  +10%.  Next  we 
select  13  stations  with  approximately  1  km  between  stations,  and  a  center  frequency  of  2 
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Hz.  We  should  get  the  same  result  as  in  the  first  case.  Instead  we  find  inconsistencies  in 
the  phase  velocities  determined  using  the  different  frequency  bands.  The  phase  velocities 
we  would  expect,  based  on  an  lASWl  crust,  are  Pn-8.1  km/sec,  P*-6.5  km/sec,  PmP-6.6 
Km/sec,  Pg-5.8  km/sec,  Sn-4.5  km/sec,  SmS-3.8  km/sec.  If  we  assume  a  crustal  model 
such  as  that  of  Kanestrom  and  Haugland  (1971),  as  described  by  VogQord  and  Langston 
(1990),  we  expect  slightly  higher  phase  velocities.  Here,  for  example,  we  should  find  Pn- 
8.2  km/sec,  PmP-7.3  km/sec,  Sn-4.8  km/sec,  SmS-4.2  lan/sec.  We  see  that  the  F-K  plots 
give  unreasonable  results  for  the  first  P  airiv^  (Pn)  at  both  frequencies,  and  questionable 
results  for  other  phases.  Note  that  the  SmS  velocity  is  getting  low  enough  that  the 
corresponding  wavenumber  is  at  the  borderline  of  teing  acceptable  for  analysis  with  these 
values  of  station  spacing,  due  to  possible  aliasing. 

RECOMMENDATIONS  AND  FUTURE  PLANS 

It  is  clear  from  the  above  examples  that  a  combination  of  raytracing  and  waveform 
modeling  can  be  an  effective  way  to  characterize  regions  with  few  events  and  fewer 
stations.  It  must  however  be  emphasized  that  waveform  modeling  without  ground  truth  is 
subject  to  misinterpretation.  One  important  factor  is  the  misidentification  of  phases.  This 
has  been  mentioned  by  Vogfjord  and  Langston  (1990).  Phases  may  be  missed  because 
they  have  low  amplitude,  or  mislabeled  b^ause  of  simplified  assumptions.  For  example  at 
short  distances  (less  than  150-200  km)  Lg  is  not  a  phase  with  a  constant  group  velocity  of 
3.5  km/sec.  One  can  instead  look  for  the  SmS  phase,  which  can  be  quite  strong  and  which 
comprises  the  beginning  of  the  Lg  wavetrain  at  greater  distances.  The  second  factor  is  the 
bias  provided  by  the  crust  and  upper  mantle.  This  can  be  due  to  anisotropy  or  to  three 
dimensional  effects.  For  example  one  factor  that  may  contribute  to  the  difficulty  in 
determining  reasonable  phase  velocities  for  the  Bl^j0  (Quarry  event  mentioned  above  is  the 
presence  of  the  Oslo  Graben  southwest  of  the  NORESS  array  (see  graben  location  in 
Vogfjord  and  Langston,  1990). 

Further  work  is  required  to  quantify  the  effects  of  geologically  realistic 
heterogeneities  on  travel  times  and  amplitudes.  The  ability  to  do  rapid  relocations  using  a 
variety  of  velocity  models  through  a  range  of  near  to  far  regional  distances  is  important, 
and  we  are  now  focusing  our  efforts  in  this  area.  We  are  basing  our  travel  time  calculations 
on  the  tau-p  method  of  Buland  and  Chapman  (1993),  for  which  computation  is  rapid.  If 
we  do  not  have  detailed  independent  velocity  information  in  an  area  it  seems  reasonable  to 
work  with  a  starting  one-dimensional  crustal  model,  and  quantify  how  variations  in  that 
model  would  affect  our  locations.  It  seems  clear  that  a  simple  global  model  like  IASP91  is 
unlikely  to  be  sufficient  in  specific  regions  for  waveform  modeling  and  reliable  locations 
within  a  500  km  range. 

Use  of  arrays  for  regional  phase  velocity  determinations  is  problematic.  Previous 
work  at  NORSAR  has  suggested  Aat  P  and  S  energy  can  be  differentiated  on  F-K  plots  at 
NORESS,  but  that  it  may  sometimes  be  difficult  to  separate  out  phases  more  precisely  from 
inferred  phase  velocities  (Mykkeltveit  and  Bungum,  1984;  Mykkeltveit  et  al.,  1990).  It 
may  be  that  calibration  of  each  array  is  necessary  as  a  function  of  distance  and  frequency. 
This  is  clearly  not  appropriate  in  the  situation  described  here  in  which  little  if  any  ground 
truth  is  presumed  available,  and  in  which  we  would  like  to  use  arrays  to  constrain  crustal 
structure  and  identify  specific  P  and  S  arrivals  within  the  context  of  that  structure. 
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Figure  2.  Linear  finite  difference  synthetics  overlain  by  travel  time  curves.  Left:  Model  A  (IASP91).  Right:  Model  B  (IASP91  with 
modified  4  layer  crust).  Travel  time  curves  from  raytracing  include  direct  rays,  refractions  and  reflections.  For  example 
Model  A  includes  Pn  (#3),  P*  (2),  Pg  (#1),  P^^P  (#7),  PmP  (#8),  Sn  (#6),  S*(#5),  Sg(#4),  S^oS  (#10)  and  SmS  (#11). 
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Figure  3.  Reflectivity  velocity  synthetics  overlain  by  travel  time  curves  for  Model  B. 
Source  is  fixed  at  0  km  depth. 
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Figure  4a.  NORESS  recording  of  BlSsj0  quarry  blast  F-K  diagrams  on  the  right  are  for  signal  in  3  second  window  shown  on  left. 

Upper  right;  frequency  band  0.9- 1.1  Hz.  Six  stations  spaced  roughly  2.  km  apart  (stations  shown  on  left).  Lower  right: 
Frequency  band  1.8-2.2  Hz.  13  stations  used  here  (A  ring,  3  from  C  ring,  entire  D  ring),  spaced  roughly  1  km  apart. 
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Figure  4b.  Continuation  of  F-K  plots  from  previous  figure.  Upper  set  of  plots:  1  Hz  center  frequency,  6  stations.  Lower  set  of 
plots.  2  Hz  center  frequency,  13  stations.  One  P  and  two  S  windows  shown.  Refer  to  Figure  4a  for  location  of 


